This document presents preliminary analyses characterizing the NDVI of certain places and census tracts throughout the Denver area.

1 NDVI for the Denver area, 2016-2021

We obtained the normalized difference vegetation index (NDVI) from the Landsat-8 satellite at a resolution of 30 meters squared at a roughly 15-day interval between 2016 and 2021 for the following area around Denver:

1.1 Define study area

We gathered this data using the rgee package, an R package that facilitates connection to the Google Earth Engine via Python. Please see these two scripts for specific details on the process for gathering and processing the NDVI data in this area:

https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/1a_get_landsat_ndvi_denver.R

https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/1b_wrangle_check_landsat_ndvi_denver.R

1.2 NDVI for one day

Here is the NDVI over the study area on a cloud-free day, July 4, 2021:

#Update 5/26/22 this may not work with mapview anymore, so use plot instead for speed and
setwd(here("data-processed"))
pal_terrain_col = rev(terrain.colors(100)) 
ndvi_den_metro_terr_5_yr = terra::rast("ndvi_den_metro_terr_5_yr.tif") 

ndvi_den_metro_terr_5_yr$`20210704_NDVI` %>% plot()

#old code:
# mv_good_date =ndvi_den_metro_terr_5_yr$`20210704_NDVI` %>% 
#   raster::raster() %>%  
#   mapview(
#     col.regions = pal_terrain_col,
#     layer.name = "NDVI, July 4, 2021")
# mv_good_date

2 Characterize NDVI of places of interest

2.1 First, determine valid (cloud-free) dates

2.1.1 Testing a few places with expected high NDVI

First, pick a few places where we would expect NDVI to be high in the summer on a cloud-free day. If it’s not high, then we can assume the image is bad (i.e., clouds in the way). These test places were determined by examining the cloud-free day (July 4, 2021) above. We chose three areas: an area east of Evergreen, Colorado; a plot in City Park; and a plot in the Indian Tree Golf Club:

(These areas are defined here: https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/1b_wrangle_check_landsat_ndvi_denver.R)

setwd(here("data-processed"))
getwd()
## [1] "/Users/michaeldgarber/Dropbox/Work/CSU-UCSD post-doc/post-doc-research/green-space-denver-project/green-space-denver/data-processed"
load("bbox_evergreen_east.RData")
load("bbox_indian_tree_golf.RData")
load("bbox_city_park.RData")
load("den_metro_bbox_custom.RData")
mv_evergreen_east = bbox_evergreen_east %>%
  mapview(col.regions = "red",  layer.name = "Evergreen East")
mv_indian_tree_golf = bbox_indian_tree_golf %>%
  mapview(col.regions = "red", layer.name = "Indian Tree Golf")
mv_city_park =bbox_city_park %>%
  mapview(col.regions = "red",  layer.name = "City Park plot")
mv_den_metro_bbox_custom = den_metro_bbox_custom %>% 
  mapview(layer.name = "Study area", alpha.regions = .2)
mv_den_metro_bbox_custom+ mv_evergreen_east + mv_indian_tree_golf + mv_city_park

2.1.2 What is NDVI over time for these test places?

Let’s look at the temporal NDVI trend for these three places, without excluding any dates. This is the mean NDVI on that day. That is, each place includes several pixels of size 30 meters squared. This is the average of the NDVI for those pixels for that place on that day.

setwd(here("data-processed"))
load("ndvi_test_places_day_wrangle.RData")
ndvi_test_places_day_wrangle %>% 
  ggplot(
    aes(
      x=date, 
      y=ndvi_mean  
    ))+
  geom_line(aes(colour = test_place_name))+
  geom_point(aes(colour = test_place_name))+
  scale_x_date(labels=date_format("%Y-%m-%d"), date_breaks = "3 months")+
  theme(axis.text.x = element_text(angle = 45, hjust=1))

Hmm, if we look closely, we can see that there are some very low NDVI values for some of these places even in the summer, which is not plausible. Those low values must indicate obstruction by clouds. For example, look at a day when NDVI was measured as very low in City Park despite it being mid-May when we would expect it to be higher. So I’m going to exclude dates like these.

setwd(here("data-processed"))
pal_viridis_trunc=viridis::viridis_pal(end=.9) #trunc for truncated
ndvi_den_metro_terr_5_yr = terra::rast("ndvi_den_metro_terr_5_yr.tif")
ndvi_den_metro_terr_5_yr$`20180517_NDVI` %>% #Use 2018-05-17
  plot()

#old code that uses mapview. decide to use plot instead. see abobve
# mv_bad_date =ndvi_den_metro_terr_5_yr$`20180517_NDVI` %>% #Use 2018-05-17
#   raster::raster() %>%  
#   mapview(
#     col.regions = pal_viridis_trunc,
#     layer.name = "NDVI, May 17, 2018")
# mv_bad_date

The specific NDVI thresholds for each thresholds are noted in this script (https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/1b_wrangle_check_landsat_ndvi_denver.R) and warrant further discussion.

The valid dates are restricted to May, June, July, and August and specifically are:

setwd(here("data-processed"))
load("lookup_date_is_valid_all.RData") #load dataset of valid dates created elsewhere
lookup_date_is_valid_all %>% 
  filter(date_is_valid_all==1) %>% 
  dplyr::select(date) %>% 
  pull()
##  [1] "2016-05-08" "2016-06-01" "2016-06-09" "2016-06-17" "2016-07-03"
##  [6] "2016-08-12" "2017-05-09" "2017-06-02" "2017-06-26" "2017-07-12"
## [11] "2017-07-28" "2017-08-13" "2017-08-21" "2017-08-29" "2018-05-25"
## [16] "2018-06-10" "2018-06-18" "2018-06-26" "2018-07-04" "2018-07-12"
## [21] "2018-08-05" "2018-08-13" "2018-08-21" "2018-08-29" "2019-05-09"
## [26] "2019-06-26" "2019-07-04" "2019-08-21" "2019-08-29" "2020-05-24"
## [31] "2020-06-25" "2020-07-27" "2020-08-12" "2020-08-28" "2021-05-25"
## [36] "2021-06-02" "2021-06-10" "2021-07-04" "2021-07-28" "2021-08-29"

2.2 Load places of interest

Load the places-of-interest dataset, which is created here: https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/0_read_denver_native_zones.R

setwd(here("data-processed"))
load("places_native_geo.RData")
load("places_native_nogeo.RData")
load("native_places_ndvi_day_nogeo.RData")

Print information about the places of interest, including approximate percent nativity where available.

places_native_nogeo %>% 
  dplyr::select(place_name, native_percent, place_type, place_id) %>% 
  kable()
place_name native_percent place_type place_id
Denver Botanic Gardens, 100% Native 100 native spectrum 1
Plains Conservation Center, 100% Native 100 native spectrum 2
Green Mountain Park, 85% Native 85 native spectrum 3
Hogback along C-470, 75% Native 75 native spectrum 4
Suburban Open Space 1, Chatfield H.S., 50% Native 50 native spectrum 5
Suburban Open Space 2, Columbine Hills Church, 30% Native 30 native spectrum 6
Denver Botanic Gardens at Chatfield, 10% Native 10 native spectrum 7
City Park, 0% Native 0 native spectrum 8
Chatfield Meadow Restoration NA high diversity 9
Chatfield Prairie Garden NA high diversity 10
City Park Greenhouses NA high diversity 11
Denver Botanic Gardens Green Roof NA high diversity 12
Kenderick Lake Xeriscape Garden NA high diversity 13

2.2.1 Map places of interest by type

Map these places by type (nativity spectrum or high-plant diversity).

2.3 Graph NDVI over time for these places

During valid dates only (above).

2.3.1 For the places with percent nativity values

The first set of polygons sent over.

native_places_ndvi_day_nogeo %>% 
  filter(date_is_valid_all==1) %>% 
  filter(place_type == "native spectrum") %>% 
  ggplot(aes(
    x=date, 
    y=ndvi_med # median
  ))+
  geom_ribbon( #ribbon around 25th and 75th percentile
    aes(ymin =ndvi_25, ymax = ndvi_75 ),
    alpha=.4
  )+
  ylab("NDVI, Median") +
  xlab("Date") +
  geom_line( size=.7 ) +#note size better than lwd
  geom_point()+
  scale_x_date(breaks = pretty_breaks())+
  scale_y_continuous(
    limits = c(0, NA),  #force axis origin to be zero
    breaks= seq(0,0.8,by = 0.1))+ 
  theme(
    axis.text.x=element_text(angle=60, hjust=1),
        panel.border = element_rect(
          colour = "gray72", size = 0.5, fill=NA))+
  facet_grid(   # facet them
    cols = vars(place_name_fac),
    labeller = label_wrap_gen() #wrap facet labels
  )

2.3.2 For the places with high plant diversity

This is the second set of polygons sent.

native_places_ndvi_day_nogeo %>% 
  filter(date_is_valid_all==1) %>% 
  filter(place_type == "high diversity") %>% 
  ggplot(aes(
    x=date, 
    y=ndvi_med # median
  ))+
  geom_ribbon( 
    aes(ymin =ndvi_25, ymax = ndvi_75 ),
    alpha=.4
  )+
  ylab("NDVI, Median") +
  xlab("Date") +
  geom_line( size=.7 ) + 
  geom_point()+
  scale_x_date(breaks = pretty_breaks())+
  scale_y_continuous(
    limits = c(0, NA),  
    breaks= seq(0,0.8,by = 0.1))+ 
  theme(
    axis.text.x=element_text(angle=60, hjust=1),
        panel.border = element_rect(
          colour = "gray72", size = 0.5, fill=NA))+
  facet_grid(   
    cols = vars(place_name_fac),
    labeller = label_wrap_gen() #wrap facet labels
  )

2.4 Graph NDVI against percent nativity

Among the first set of polygons with values for percent nativity.

native_places_ndvi_day_nogeo %>% 
  filter(date_is_valid_all==1) %>% 
  filter(place_type == "native spectrum") %>% 
  ggplot(aes(
    x=native_percent, 
    y=ndvi_med # median
  ))+
  geom_point(
    aes(colour=place_name_fac), size = 1.5, 
    alpha=.5 #varying alpha to illustrate density
    ) +
  xlab("Percent native") +
  ylab("NDVI, median") +
  scale_y_continuous(
    limits = c(0, NA),  #force axis origin to be zero
    breaks= seq(0,0.8,by = 0.1))+ 
  scale_color_hue(
    name = "Place name"
  )

2.5 Summary of NDVI for these places

setwd(here("data-processed"))
load("places_native_nogeo.RData")
load("native_places_ndvi_day_nogeo.RData")
native_places_ndvi_day_nogeo %>% 
  filter(date_is_valid_all==1) %>% 
  group_by(place_type, place_name_fac) %>% 
  summarise(
    ndvi_mean = mean(ndvi_mean_wt, na.rm=TRUE), # mean of means
    ndvi_25th = quantile(ndvi_med, probs =c(0.25), na.rm=TRUE), #percentile of medians
    ndvi_med = median(ndvi_med, na.rm=TRUE), #percentile of medians
    ndvi_75th = quantile(ndvi_med, probs =c(0.75), na.rm=TRUE) ) %>% #percentile of medians
  ungroup() %>% 
  left_join(places_native_nogeo, by = c("place_name_fac", "place_type")) %>% 
  dplyr::select(
    starts_with("place_type"),
    starts_with("place_name_fac"),
    starts_with("area_ft2"),
    starts_with("area_mi"),
    starts_with("ndvi_mea"),
    starts_with("ndvi_25"),
    starts_with("ndvi_med"),
    starts_with("ndvi_75th"),
  ) %>% 
  as_tibble() %>% 
  knitr::kable(
    booktabs = TRUE,
    col.names = c("Place type", "Place name", "Area (sq. ft.)", "Area (sq. mi.)", "NDVI, mean", "NDVI, 25th-ile", "NDVI, 50th-ile", "NDVI, 75th-ile"),
  digits = c(0, 0, 0, 3, 2, 2, 2, 2, 2)
  )
Place type Place name Area (sq. ft.) Area (sq. mi.) NDVI, mean NDVI, 25th-ile NDVI, 50th-ile NDVI, 75th-ile
high diversity Chatfield Meadow Restoration 4391 0.000 0.48 0.47 0.50 0.50
high diversity Chatfield Prairie Garden 2922 0.000 0.37 0.35 0.37 0.37
high diversity City Park Greenhouses 9596 0.000 0.31 0.27 0.28 0.28
high diversity Denver Botanic Gardens Green Roof 12882 0.000 0.34 0.34 0.36 0.36
high diversity Kenderick Lake Xeriscape Garden 35418 0.001 0.40 0.36 0.37 0.37
native spectrum Denver Botanic Gardens, 100% Native 21585 0.001 0.54 0.52 0.55 0.55
native spectrum Plains Conservation Center, 100% Native 793905 0.028 0.16 0.23 0.26 0.26
native spectrum Green Mountain Park, 85% Native 83536011 2.996 0.42 0.36 0.40 0.40
native spectrum Hogback along C-470, 75% Native 1293880 0.046 0.36 0.32 0.35 0.35
native spectrum Suburban Open Space 1, Chatfield H.S., 50% Native 161894 0.006 0.40 0.37 0.41 0.41
native spectrum Suburban Open Space 2, Columbine Hills Church, 30% Native 198596 0.007 0.36 0.26 0.31 0.31
native spectrum Denver Botanic Gardens at Chatfield, 10% Native 97135 0.003 0.53 0.40 0.54 0.54
native spectrum City Park, 0% Native 2664873 0.096 0.67 0.67 0.70 0.70

3 NDVI of census tracts

In this section, we summarize NDVI for census tracts in Denver County. Before characterizing NDVI of each census tract, we removed bodies of water. NDVI values of water approach -1, and we did not want to “penalize” a census tract for having large amounts of water.

Here is a map of the median NDVI of the census tracts in the area on 2021-07-04:

Here is a histogram of the distribution of median NDVI for these census tracts on June 10, 2021:

den_co_tract_ndvi_day_geo %>% 
  filter(date == "2021-06-10") %>% 
  ggplot(aes(ndvi_med))+
  geom_histogram()

4 NDVI of public green spaces

In this section, we examine the NDVI levels of public green space in this area, including city parks and nature reserves. Please see this code for more details on how we gathered the polygons for the greens spaces: https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/0_load_denver_osm_parks_water.R

Here are the parks we downloaded:

setwd(here("data-processed"))
load("den_jeff_co_green_space_public.RData")
den_jeff_co_green_space_public %>% 
  mapview(zcol = "osm_value", layer.name = "Green space type")

Like with census tracts, we removed bodies of water from these green spaces before measuring NDVI. Here is the NDVI on June 10, 2021:

setwd(here("data-processed"))
load("den_metro_green_space_ndvi_day_geo.RData")
load("den_jeff_co_geo.RData")
pal_viridis_trunc=viridis::viridis_pal(end=.9) #trunc for truncated
mv_den_jeff_co_geo = den_jeff_co_geo %>%  
  mapview(
    layer.name = "County name",
    color = c("red", "orange"),
    col.regions = c("red", "orange"),
    lwd=2,
    zcol = "county_name_short", alpha.regions = 0)
mv_den_metro_green_space_ndvi_day_geo = den_metro_green_space_ndvi_day_geo %>%
  dplyr::select( #limit some of the variables before mapview
    starts_with("osm"),
    starts_with("ndvi"),
    starts_with("date"),
    starts_with("county")) %>% 
  filter(date == "2021-06-10") %>% 
  mapview(
    layer.name = "NDVI, Median",
    col.regions = pal_viridis_trunc,
    zcol = "ndvi_med" 
  )
mv_den_metro_green_space_ndvi_day_geo+mv_den_jeff_co_geo

Most parks are in the city, and most protected areas are in the mountains. Summarize NDVI by type of green space:

setwd(here("data-processed"))
library(knitr)
load("den_metro_green_space_ndvi_day_nogeo.RData")
den_metro_green_space_ndvi_day_nogeo %>% 
  filter(date_is_valid_all==1) %>% 
  group_by(osm_value) %>% 
  summarise(
    ndvi_25th = quantile(ndvi_med, probs =c(0.25), na.rm=TRUE),
    ndvi_med = median(ndvi_med, na.rm=TRUE),
    ndvi_75th = quantile(ndvi_med, probs =c(0.75), na.rm=TRUE),
    ) %>% 
  ungroup() %>% 
  knitr::kable()
osm_value ndvi_25th ndvi_med ndvi_75th
nature_reserve 0.3313458 0.4128266 0.4128266
park 0.3722610 0.4735859 0.4735859
protected_area 0.3428352 0.4256270 0.4256270


Copyright © 2022 Michael D. Garber